#include <string>
#include <time.h>
#include <numeric>
#include <stdlib.h>

#include "BaseAlgorithm.h"
#include "DataTypeDefine.h"

// itk include
#include "itkMatrix.h"
#include "itkSample.h"
#include "itkAddImageFilter.h"
#include "itkMultiplyImageFilter.h"
#include "itkThresholdImageFilter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "itkMinimumMaximumImageCalculator.h"
#include "itkAdditiveGaussianNoiseImageFilter.h"
#include "itkDanielssonDistanceMapImageFilter.h"
#include "itkSignedMaurerDistanceMapImageFilter.h"
#include "itkSmoothingRecursiveGaussianImageFilter.h"

#include "itkSimilarity3DTransform.h"
#include "itkEuler3DTransform.h"
#include "itkAffineTransform.h"
#include "itkScaleTransform.h"

// vtk include
//#include <vtkMatrix4x4.h>
//#include <vtkSmartPointer.h>
//#include <vtkTransform.h>

#define MAX_NUMBER 28
#define SAMPLE_PER_VERTEBRA 2  // 每节椎体的采样数量
#define BOUNDING_SIZE 128
#define ROTATE_ANGLE 12
#define SHRINK_SIZE 2

#define boundingSizeX 96
#define boundingSizeY 96
#define boundingSizeZ 80

#define Random(x) (rand()%x)  // [0, x)
#define Random2(x) (Random(2 * (x) + 1) - (x))  // [-x, x]

template<class S, class T> inline T min(const S &a, const T &b) { if (a<b) return (T)a; else return b; }
template<class S, class T> inline T max(const S &a, const T &b) { if (a>b) return (T)a; else return b; }

/**
* @brief 获取图像的最大最小值
* @param image    输入的图像
* @param minValue 图像的最小值
* @param maxValue 图像的最大值
*
*/
void minMaxFilter(LabelMapPointer image, int &minValue, int &maxValue)
{
	typedef itk::ImageRegionIteratorWithIndex<LabelMapType> imageIteratorType;
	imageIteratorType imageIterator(image, image->GetBufferedRegion());

	minValue = 25;
	maxValue = 0;

	while (!imageIterator.IsAtEnd())
	{
		const int value = imageIterator.Get();
		++imageIterator;

		if (value == 0)
		{
			continue;
		}

		minValue = std::min(value, minValue);
		maxValue = std::max(value, maxValue);
	}

	if (minValue == 25)
	{
		minValue = 0;
	}
}

/**
* @brief 对输入图像做乘法
* @param image    输入的图像
* @param coeff    乘数
*
* @return 返回输出图像
*/
LabelMapPointer multiplyLabel(LabelMapPointer image, int coeff)
{
	typedef itk::MultiplyImageFilter<LabelMapType, LabelMapType> MultiplyImageFilterType;
	MultiplyImageFilterType::Pointer multiplyImageFilter = MultiplyImageFilterType::New();
	multiplyImageFilter->SetInput(image);
	multiplyImageFilter->SetConstant(coeff);
	multiplyImageFilter->Update();
	return multiplyImageFilter->GetOutput();
}

/**
* @brief 获取图像可裁剪范围
* @param image 输入的图像
* @param range 图像可裁剪范围
* @return 返回获取结果，false表示图像尺寸过小，无法裁剪
*/
bool getClipRange(CTImagePointer image, int* clipRange)
{
	CTImageType::SizeType dimensions = image->GetLargestPossibleRegion().GetSize();

	if (dimensions[0] < boundingSizeX || dimensions[1] < boundingSizeY || dimensions[2] < boundingSizeZ)
	{
		return false;
	}

	clipRange[0] = dimensions[0] - boundingSizeX;
	clipRange[1] = dimensions[1] - boundingSizeY;
	clipRange[2] = dimensions[2] - boundingSizeZ;

	return true;
}

/**
* @brief 对输入图像进行padding，具体方向请参考itk坐标系
* @param image    输入的图像
* @param I        I方向padding的大小
* @param S        S方向padding的大小
* @param L		  L方向padding的大小
* @param R		  R方向padding的大小
* @param P        P方向padding的大小
* @param A        A方向padding的大小
*
* @return padding后的图像
*/
template<typename Tin> 
typename Tin::Pointer extendImage(typename Tin::Pointer image, int I, int S, int L, int R, int P, int A)
{
	typename Tin::PointType origin = image->GetOrigin();
	typename Tin::DirectionType direction = image->GetDirection();
	typename Tin::SpacingType spacing = image->GetSpacing();
	typename Tin::SizeType inputSize = image->GetLargestPossibleRegion().GetSize();

	typename Tin::RegionType region;
	typename Tin::IndexType start;
	start[0] = 0;
	start[1] = 0;
	start[2] = 0;

	typename Tin::SizeType size;
	size[0] = inputSize[0] + L + R;
	size[1] = inputSize[1] + P + A;
	size[2] = inputSize[2] + I + S;

	region.SetSize(size);
	region.SetIndex(start);

	// 分配空间
	typename Tin::Pointer extended = Tin::New();
	extended->SetRegions(region);
	extended->Allocate();
	extended->SetSpacing(spacing);
	extended->SetOrigin(origin);
	extended->SetDirection(direction);

	for (unsigned int row = 0; row < size[0]; row++)
	{
		for (unsigned int col = 0; col < size[1]; col++)
		{
			for (unsigned int h = 0; h < size[2]; h++)
			{
				typename Tin::IndexType tgtIndex;
				typename Tin::IndexType srcIndex;
				tgtIndex[0] = row;
				tgtIndex[1] = col;
				tgtIndex[2] = h;
				if ((h   > I) && (h   < size[2] - S) &&
					(row > L) && (row < size[0] - R) &&
					(col > P) && (col < size[1] - A))
				{
					srcIndex[0] = row - L;
					srcIndex[1] = col - P;
					srcIndex[2] = h - I;
					extended->SetPixel(tgtIndex, image->GetPixel(srcIndex));
				}
				else
				{
					extended->SetPixel(tgtIndex, 0);
				}
			}
		}
	}

	return extended;

}

/**
* @brief 对输入图像中一定范围内的标记值分别进行距离变换
* @param image       输入的图像
* @param startLabel  起始标记值
* @param endLabel    终止标记值
*
* @return 返回距离图的数组
*/
std::vector<ImagePointer> distanceTransform(LabelMapPointer image, int startLabel, int endLabel)
{
	typedef itk::SignedMaurerDistanceMapImageFilter<LabelMapType, ImageType> SignedMaurerDistanceMapImageFilterType;
	
	std::vector<ImagePointer> distanceMaps(endLabel - startLabel + 1, ImageType::New());

	for (int label = startLabel; label <= endLabel; ++label)
	{
		// std::cout << label << std::endl;
		LabelMapPointer oneLabel = BaseAlgorithm::ThresholdImage<LabelMapType>(image, label * 100, label * 100, 1);

		SignedMaurerDistanceMapImageFilterType::Pointer distanceMapImageFilter = SignedMaurerDistanceMapImageFilterType::New();
		distanceMapImageFilter->SetInput(oneLabel);
		distanceMapImageFilter->Update();
		distanceMaps[label - startLabel] = distanceMapImageFilter->GetOutput();
	}

	return distanceMaps;
}

/**
* @brief 对输入图像进行距离变换
* @param image       输入的图像
*
* @return 返回距离图
*/
ImagePointer distanceTransform(LabelMapPointer image)
{
	typedef itk::SignedMaurerDistanceMapImageFilter<LabelMapType, ImageType> SignedMaurerDistanceMapImageFilterType;

	SignedMaurerDistanceMapImageFilterType::Pointer distanceMapImageFilter = SignedMaurerDistanceMapImageFilterType::New();
	distanceMapImageFilter->SetInput(image);
	distanceMapImageFilter->Update();

	return distanceMapImageFilter->GetOutput();
}

/**
* @brief 对输入图像进行预处理
* @param image       输入的图像
* @oaram minLabel    通过引用返回标记最小值
* @param maxLabel    通过引用返回标记最大值
*
* @return 返回处理结果，每个结果由大小为10的数组bounds储存，形式如下
*     bounds[0]:             每节椎体的体积
*     bounds[1]~bounds[6]:   每节椎体的包围盒大小
*     bounds[7]~bounds[9]:   每节椎体的中心位置
*/
std::vector<std::vector<int>> prepareData(LabelMapPointer labelmap, int& minLabel, int& maxLabel)
{
	std::vector<std::vector<int>> bounds(MAX_NUMBER, std::vector<int>(10, 0));

	for (int id = 0; id < MAX_NUMBER; id++)
	{
		bounds[id][1] = 10000;
		bounds[id][3] = 10000;
		bounds[id][5] = 10000;
	}

	LabelMapType::IndexType index;
	int value;

	minLabel = MAX_NUMBER;
	maxLabel = 0;

	typedef itk::ImageRegionIteratorWithIndex<LabelMapType> imageIteratorType;
	imageIteratorType iterator(labelmap, labelmap->GetLargestPossibleRegion());
	while (!iterator.IsAtEnd())
	{
		value = iterator.Get();
		index = iterator.GetIndex();

		if (value > 0)
		{
			bounds[value][7] += index[0];
			bounds[value][8] += index[1];
			bounds[value][9] += index[2];

			minLabel = std::min(minLabel, value);
			maxLabel = std::max(maxLabel, value);

			bounds[value][0]++;
			bounds[value][1] = min(index[0], bounds[value][1]);
			bounds[value][2] = max(index[0], bounds[value][2]);
			bounds[value][3] = min(index[1], bounds[value][3]);
			bounds[value][4] = max(index[1], bounds[value][4]);
			bounds[value][5] = min(index[2], bounds[value][5]);
			bounds[value][6] = max(index[2], bounds[value][6]);
		}
		++iterator;
	}

	for (int i = 0; i < MAX_NUMBER; i++)
	{
		if (bounds[i][0] > 0)
		{
			bounds[i][7] /= bounds[i][0];
			bounds[i][8] /= bounds[i][0];
			bounds[i][9] /= bounds[i][0];
		}
	}

	return bounds;
}

/**
* @brief 对输入图像增加高斯噪声
* @param image       输入的图像
*
* @return 返回处理结果
*/
ImagePointer gaussNoise(ImagePointer image)
{
	typedef itk::AdditiveGaussianNoiseImageFilter<ImageType, ImageType> AddNoiseFilterType;
	AddNoiseFilterType::Pointer addNoiseFilter = AddNoiseFilterType::New();

	addNoiseFilter->SetInput(image);
	addNoiseFilter->SetMean(20 + Random(5));
	addNoiseFilter->SetStandardDeviation(20 + Random(5));
	addNoiseFilter->Update();

	return addNoiseFilter->GetOutput();
}

/**
* @brief 对输入图像增加高斯平滑
* @param image       输入的图像
*
* @return 返回处理结果
*/
ImagePointer gaussSmoothing(ImagePointer image)
{
	typedef itk::SmoothingRecursiveGaussianImageFilter<ImageType, ImageType> SmoothingFilterType;
	SmoothingFilterType::Pointer smoothingFilter = SmoothingFilterType::New();

	smoothingFilter->SetInput(image);
	smoothingFilter->SetSigma(0.25);
	smoothingFilter->Update();

	return smoothingFilter->GetOutput();
}

/**
* @brief 对输入图像进行翻转
* @param image       输入的图像
* @param axis        翻转的维度
*
* @return 返回处理结果
*/
CTImagePointer flipImage(CTImagePointer image, int axis)
{
	typedef itk::FlipImageFilter<CTImageType> FlipImageFilterType;
	FlipImageFilterType::Pointer flipImageFilter = FlipImageFilterType::New();

	FlipImageFilterType::FlipAxesArrayType flipAxes;
	flipAxes[0] = false;
	flipAxes[1] = false;
	flipAxes[2] = false;
	flipAxes[axis] = true;

	flipImageFilter->SetInput(image);
	flipImageFilter->SetFlipAxes(flipAxes);
	return flipImageFilter->GetOutput();
}

/**
* @brief 对输入图像的像素增加一个定值
* @param image       输入的图像
* @param add         增加的值
*
* @return 返回处理结果
*/
CTImagePointer addImage(CTImagePointer image, int add)
{
	typedef itk::AddImageFilter<CTImageType> AddImageFilterType;
	AddImageFilterType::Pointer addImageFilter = AddImageFilterType::New();

	addImageFilter->SetInput1(image);
	addImageFilter->SetConstant2(add);
	addImageFilter->Update();

	return addImageFilter->GetOutput();
}

/**
* @brief 对输入图像进行旋转
* @param image             输入的图像
* @param angelX	    	   X轴的旋转值
* @param angelY            Y轴的旋转值
* @param angelZ            Z轴的旋转值
* @param interpolateType   插值的形式，有以下几种选择：
*     0 线性插值
*     1 最近邻插值
*
* @return 返回处理结果
*/
template<typename Tin>
typename Tin::Pointer transformImage(typename Tin::Pointer image, double angleX, double angleY, double angleZ, int interpolateType = 0)
{
	typename Tin::PointType origin = image->GetOrigin();
	typename Tin::DirectionType direction = image->GetDirection();
	typename Tin::SpacingType spacing = image->GetSpacing();
	typename Tin::SizeType inputSize = image->GetLargestPossibleRegion().GetSize();
	typename Tin::IndexType corner = image->GetLargestPossibleRegion().GetIndex();

	// 获取实际的原点（corner不一定等于origin）
	origin[0] += direction.GetVnlMatrix().get(0, 0) * spacing[0] * corner[0];
	origin[1] += direction.GetVnlMatrix().get(1, 1) * spacing[1] * corner[1];
	origin[2] += direction.GetVnlMatrix().get(2, 2) * spacing[2] * corner[2];

	// 获取实际的旋转中心
	double centerX = origin[0] + direction.GetVnlMatrix().get(0, 0) * spacing[0] * boundingSizeX / 2;
	double centerY = origin[1] + direction.GetVnlMatrix().get(1, 1) * spacing[1] * boundingSizeY / 2;
	double centerZ = origin[2] + direction.GetVnlMatrix().get(2, 2) * spacing[2] * boundingSizeZ / 2;

	typedef itk::Euler3DTransform<double> eularType;
	eularType::Pointer rotation = eularType::New();
	eularType::InputPointType center;
	center[0] = centerX;
	center[1] = centerY;
	center[2] = centerZ;
	rotation->SetRotation(angleX, angleY, angleZ);
	rotation->SetCenter(center);

	typedef	itk::ResampleImageFilter<Tin, Tin> resampleFilterType;
	typename resampleFilterType::Pointer rotateResampleFilter = resampleFilterType::New();
	typename resampleFilterType::Pointer scaleResampleFilter = resampleFilterType::New();

	if (interpolateType == 0)
	{
		typedef itk::LinearInterpolateImageFunction<Tin, double> LinearInterpolatorType;
		typename LinearInterpolatorType::Pointer linearInterpolator = LinearInterpolatorType::New();
		rotateResampleFilter->SetInterpolator(linearInterpolator);
		scaleResampleFilter->SetInterpolator(linearInterpolator);
	}
	else
	{
		typedef itk::NearestNeighborInterpolateImageFunction<Tin, double> NNInterpolatorType;
		typename NNInterpolatorType::Pointer nearersInterpolator = NNInterpolatorType::New();
		rotateResampleFilter->SetInterpolator(nearersInterpolator);
		scaleResampleFilter->SetInterpolator(nearersInterpolator);
	}

	rotateResampleFilter->SetInput(image);
	rotateResampleFilter->SetTransform(rotation);
	rotateResampleFilter->SetSize(inputSize);
	rotateResampleFilter->SetOutputSpacing(spacing);
	rotateResampleFilter->SetOutputOrigin(origin);
	rotateResampleFilter->SetOutputDirection(direction);
	rotateResampleFilter->Update();

	typedef itk::ScaleTransform<double> scaleType;
	scaleType::Pointer scale = scaleType::New();
	scale->SetCenter(center);
	scale->SetScale(0.92);

	scaleResampleFilter->SetInput(rotateResampleFilter->GetOutput());
	scaleResampleFilter->SetTransform(scale);
	scaleResampleFilter->SetSize(inputSize);
	scaleResampleFilter->SetOutputSpacing(spacing);
	scaleResampleFilter->SetOutputOrigin(origin);
	scaleResampleFilter->SetOutputDirection(direction);
	scaleResampleFilter->Update();

	return scaleResampleFilter->GetOutput();

}

/**
* @brief 对输入图像的像素增加一个定值
* @param image       原始图像
* @param memory      记忆图像
* @param labelmap    标记图像
* @param distance    距离图像
* @param count       当前保存的数量，会进行自增
*
*/
void writeResult(ImagePointer image, LabelMapPointer memory, LabelMapPointer labelmap, ImagePointer distance, int &count)
{
	std::string basePath = "D:/data/dataset2/";
	std::string outputLabelMapName = basePath + "labelmap/labelmap_" + std::to_string(count) + ".nii.gz";
	std::string outputImageName    = basePath + "image/image_" + std::to_string(count) + ".nii.gz";
	std::string outputMemoryName   = basePath + "memory/memory_" + std::to_string(count) + ".nii.gz";
	std::string outputDistanceName = basePath + "distance/distance_" + std::to_string(count) + ".nii.gz";

	BaseAlgorithm::SaveImageFile<LabelMapType>(labelmap, outputLabelMapName);
	BaseAlgorithm::SaveImageFile<LabelMapType>(memory, outputMemoryName);
	BaseAlgorithm::SaveImageFile<ImageType>(distance, outputDistanceName);
	BaseAlgorithm::SaveImageFile<ImageType>(image, outputImageName);

	++count;
}

template<typename Tin>
void generateAndSave(typename Tin::Pointer image, LabelMapPointer label, int i, ImageType::SizeType lowerBound, int &count, bool isCropped = true)
{
	/*************************************************
	根据ROI区域裁剪图像或输入裁剪后的图像，翻转，增加噪声/平滑，旋转，并保存
	要修改i为当前值，而不是i+1为当前值
	*************************************************/

	ImagePointer croppedImage = ImageType::New();
	LabelMapPointer croppedLabelMap = LabelMapType::New();

	if (isCropped)
	{
		croppedImage = BaseAlgorithm::CastImage<Tin, ImageType>(image);
		croppedLabelMap = label;
	}
	else
	{
		typename Tin::SizeType dimensions = image->GetLargestPossibleRegion().GetSize();
		ImageType::SizeType upperBound;
		upperBound[0] = -int(lowerBound[0]) + dimensions[0] - boundingSizeX;
		upperBound[1] = -int(lowerBound[1]) + dimensions[1] - boundingSizeY;
		upperBound[2] = -int(lowerBound[2]) + dimensions[2] - boundingSizeZ;

		croppedImage = BaseAlgorithm::CropImage<Tin, ImageType>(image, lowerBound, upperBound);
		croppedLabelMap = BaseAlgorithm::CropImage<LabelMapType, LabelMapType>(label, lowerBound, upperBound);
	}
	// generate memory
	LabelMapPointer finalLabelMap = BaseAlgorithm::ThresholdImage<LabelMapType>(croppedLabelMap, i, i, 1);
	LabelMapPointer finalMemory;
	if (i == 1) // i == 1时要区别处理，直接创建空白图像(目前就通过threshold创建)
	{
		finalMemory = BaseAlgorithm::ThresholdImage<LabelMapType>(croppedLabelMap, 0, 26, 0, 1);
	}
	else
	{
		finalMemory = BaseAlgorithm::ThresholdImage<LabelMapType>(croppedLabelMap, 1, i - 1, 1, 0);
	}
	
	ImagePointer finalDistance = distanceTransform(finalLabelMap);
	writeResult(croppedImage, finalMemory, finalLabelMap, finalDistance, count);

	// gauss noise or gauss smoothing
	ImagePointer enhanced;
	if (Random(2))
	{
		enhanced = gaussSmoothing(croppedImage);
	}
	else
	{
		enhanced = gaussNoise(croppedImage);
	}
	writeResult(enhanced, finalMemory, finalLabelMap, finalDistance, count);

	// flip image
	int flipAxis = Random(2);
	ImagePointer flippedImage = BaseAlgorithm::FlipImage<ImageType>(croppedImage, flipAxis);
	LabelMapPointer flippedLabelMap = BaseAlgorithm::FlipImage<LabelMapType>(finalLabelMap, flipAxis);
	LabelMapPointer flippedMemory = BaseAlgorithm::FlipImage<LabelMapType>(finalMemory, flipAxis);
	finalDistance = BaseAlgorithm::FlipImage<ImageType>(finalDistance, flipAxis);
	writeResult(flippedImage, flippedMemory, flippedLabelMap, finalDistance, count);

	int useRotate = Random(2);
	if (useRotate == 1)
	{
		//std::cout << "Rotate : " << count << std::endl;
		double angleX = double(Random2(3)) / 57.3;
		double angleY = double(Random2(5) ) / 57.3;
		double angleZ = double(Random2(7)) / 57.3;
		ImagePointer rotatedImage = transformImage<ImageType>(croppedImage, angleX, angleY, angleZ);

		ImagePointer tempLabelMap = BaseAlgorithm::CastImage<LabelMapType, ImageType>(finalLabelMap);
		ImagePointer tempMemory = BaseAlgorithm::CastImage<LabelMapType, ImageType>(finalMemory);
		ImagePointer tempLabelMap2 = transformImage<ImageType>(tempLabelMap, angleX, angleY, angleZ);
		ImagePointer tempMemory2 = transformImage<ImageType>(tempMemory, angleX, angleY, angleZ);

		LabelMapPointer finalLabelMap2 = BaseAlgorithm::ThresholdImage<ImageType>(tempLabelMap2, 0.5, 1.1, 1);
		LabelMapPointer finalMemory2 = BaseAlgorithm::ThresholdImage<ImageType>(tempMemory2, 0.5, 1.1, 1);
		finalDistance = distanceTransform(finalLabelMap);

		//finalLabelMap = transformImage<LabelMapType>(finalLabelMap, angleX, angleY, angleZ, 1);
		//finalMemory = transformImage<LabelMapType>(finalMemory, angleX, angleY, angleZ, 1);
		//finalDistance = distanceTransform(finalLabelMap);
		//writeResult(rotatedImage, finalMemory, finalLabelMap, finalDistance, count);
		writeResult(rotatedImage, finalMemory2, finalLabelMap2, finalDistance, count);
	}
}

int main()
{
	int count = 720; // 保存的数量，0~719留给空白数据，另外生成

	std::string prefixName = "D:/data/trainingData/case";
	std::string suffixImageName = ".nii.gz";
	std::string suffixLabelName = "_label.nii.gz";

	int range[3] = { 0 };
	int roiCenter[3] = { 0 };
	int roiCorner[3] = { 0 };

	ImageType::SizeType lowerBound; // top left down corner of ROI
	ImageType::SizeType upperBound; // low right up corner of ROI

	CTImageType::SpacingType spacing;  // 下采样的体素大小
	spacing[0] = spacing[1] = spacing[2] = 1.25f;
	
	for (int id = 31; id > -1; --id)
	{
		printf("///////////////////// Processing image %d /////////////////////\n", id + 1);
		srand((int)time(0) * (id + 1));

		std::string imageFileName = prefixName + std::to_string(id + 1) + suffixImageName;
		std::string labelMapFileName = prefixName + std::to_string(id + 1) + suffixLabelName;

		int padding = id < 10 ? 12 : 5;

		/* 加载图像 */

		CTImagePointer inputImage = BaseAlgorithm::LoadImageFile<CTImageType>(imageFileName);
		if (id > 9)
		{
			// case1~10像素值范围是0~4071,所有图像像素值范围都统一到此
			inputImage = addImage(inputImage, 1024);
		}
		CTImagePointer resampledImage = extendImage<CTImageType>(BaseAlgorithm::ResampleImage<CTImageType>(inputImage, spacing), 
			0, 80, padding, padding, padding, padding);
		CTImageType::SizeType dimensions = resampledImage->GetLargestPossibleRegion().GetSize();

		/* 加载labelmap */

		LabelMapPointer labelMap = BaseAlgorithm::LoadImageFile<LabelMapType>(labelMapFileName);
		LabelMapPointer resampledLabelMap = extendImage<LabelMapType>(BaseAlgorithm::ResampleImage<LabelMapType>(labelMap, spacing), 
			0, 80, padding, padding, padding, padding);

		BaseAlgorithm::SaveImageFile<LabelMapType>(resampledLabelMap, "td1.nii.gz");
		BaseAlgorithm::SaveImageFile<CTImageType>(resampledImage, "td2.nii.gz");

		//system("pause");

		// 如果图像过小，直接跳过
		if (!getClipRange(resampledImage, range))
		{
			continue;
		}

		int minLabel, maxLabel;
		std::vector<std::vector<int>> bounds = prepareData(resampledLabelMap, minLabel, maxLabel);
		printf("minLabel: %d\nmaxLabel: %d\n", minLabel, maxLabel);

		//for (int i = 0; i >= minLabel; --i)
		for (int i = maxLabel; i >= minLabel; --i)
		{
			int sampled = 0; // 已经采样的数量

			// 如果是颈椎或者case1~10的最下面一节，少采样一次
			if ((i > 18) || ( i < minLabel + 1 && id > 9))
			{
				sampled -= 1;
			}

			// 如果是图像最下面一节，少采样一次
			if (i == minLabel)
			{
				sampled -= 1;
			}

			// 调试
			printf("%d %d %d %d %d %d\n", bounds[i][1], bounds[i][2], bounds[i][3], bounds[i][4], bounds[i][5], bounds[i][6]);

			while (sampled < SAMPLE_PER_VERTEBRA)
			{
				int r = int((bounds[i][6] - bounds[i][5]) / 3);
				roiCorner[0] = Random(range[0] + 1);
				roiCorner[1] = Random(range[1] + 1);
				roiCorner[2] = Random2(r / 2) + int(bounds[i][9]) - 2;
				roiCorner[2] = max(roiCorner[2] - (int)(boundingSizeZ / 2), 0);

				// 超出范围则跳过
				if (roiCorner[0] > bounds[i][1] + SHRINK_SIZE || (roiCorner[0] + boundingSizeX) < bounds[i][2] - SHRINK_SIZE ||
					roiCorner[1] > bounds[i][3] + SHRINK_SIZE || (roiCorner[1] + boundingSizeY) < bounds[i][4] - SHRINK_SIZE ||
					(roiCorner[2] + boundingSizeZ) < bounds[i][6] - SHRINK_SIZE)
				{
					printf("False: %d %d %d\n", roiCorner[0], roiCorner[1], roiCorner[2]);
					continue;
				}

				printf("True: %d %d %d\n", roiCorner[0], roiCorner[1], roiCorner[2]);

				lowerBound[0] = roiCorner[0];
				upperBound[0] = -roiCorner[0] + dimensions[0] - boundingSizeX;
				lowerBound[1] = roiCorner[1];
				upperBound[1] = -roiCorner[1] + dimensions[1] - boundingSizeY;
				lowerBound[2] = roiCorner[2];
				upperBound[2] = -roiCorner[2] + dimensions[2] - boundingSizeZ;

				ImagePointer croppedImage = BaseAlgorithm::CropImage<CTImageType, ImageType>(resampledImage, lowerBound, upperBound);
				LabelMapPointer croppedLabelMap = BaseAlgorithm::CropImage<LabelMapType, LabelMapType>(resampledLabelMap, lowerBound, upperBound);

				int maxLabelPatch{ 0 }, minLabelPatch{ 0 };
				minMaxFilter(croppedLabelMap, minLabelPatch, maxLabelPatch);

				// patch包含图像最下面的脊椎
				if (minLabelPatch == i)
				{
					maxLabelPatch = minLabelPatch + 1;
				}
				// patch包含图像最上面的脊椎
				else if (maxLabelPatch == i)
				{
					minLabelPatch = maxLabelPatch - 1;
				}
				else if (maxLabelPatch == 0)
				{
					continue;
				}
				else
				{
					minLabelPatch = i - 1;
					maxLabelPatch = i + 1;
				}

				std::cout << "flag1 " << maxLabelPatch << " " << minLabelPatch << std::endl;
				for (int j = maxLabelPatch; j >= minLabelPatch; --j)
				{
					generateAndSave<ImageType>(croppedImage, croppedLabelMap, j, lowerBound, count, true);
				}
				sampled++;
			}

			sampled = 0;
			while (sampled < SAMPLE_PER_VERTEBRA - 1 && i < 18 && i < maxLabel) 
			{
	
				roiCorner[0] = Random(range[0] + 1);
				roiCorner[1] = Random(range[1] + 1);
				roiCorner[2] = i == minLabel ? bounds[i][5] : int(bounds[i][9]) - 3;
				roiCorner[2] += Random2(5);
				roiCorner[2] = std::max(roiCenter[2], 0);

				lowerBound[0] = roiCorner[0];
				upperBound[0] = -roiCorner[0] + dimensions[0] - boundingSizeX;
				lowerBound[1] = roiCorner[1];
				upperBound[1] = -roiCorner[1] + dimensions[1] - boundingSizeY;
				lowerBound[2] = roiCorner[2];
				upperBound[2] = -roiCorner[2] + dimensions[2] - boundingSizeZ;
	
				if (roiCorner[0] > bounds[i][1] + SHRINK_SIZE * 4 || (roiCorner[0] + boundingSizeX) < bounds[i][2] - SHRINK_SIZE * 4 ||
					roiCorner[1] > bounds[i][3] + SHRINK_SIZE * 4 || (roiCorner[1] + boundingSizeY) < bounds[i][4] - SHRINK_SIZE * 4)
				{
					continue;
				}

				ImagePointer croppedImage = BaseAlgorithm::CropImage<CTImageType, ImageType>(resampledImage, lowerBound, upperBound);
				LabelMapPointer croppedLabelMap = BaseAlgorithm::CropImage<LabelMapType, LabelMapType>(resampledLabelMap, lowerBound, upperBound);

				int bias = i > 9 ? 2 : 1;
				for (int add = 0; add <= bias; ++add)
				{
					std::cout << "test 1 : " << count + 1 << std::endl;
					generateAndSave<ImageType>(croppedImage, croppedLabelMap, i + add , lowerBound, count, true);
				}
				sampled++;
			}

			// 12/2 改到这里
			sampled = 0;
			if (!(id < 10 || id == 16 || id == 27 || id == 30 || id == 29) || i > minLabel)
			{
				sampled = SAMPLE_PER_VERTEBRA; // continue
			}

			std::cout << "Phase 3: " << std::endl;

			while (sampled < SAMPLE_PER_VERTEBRA)
			{
				roiCorner[0] = Random(range[0] + 1);
				roiCorner[1] = Random(range[1] + 1);
				roiCorner[2] = int(bounds[i][9]) - (int)(boundingSizeZ / 2) - Random(20);
				roiCorner[2] = std::max(roiCorner[2], 0);

				lowerBound[0] = roiCorner[0];
				upperBound[0] = -roiCorner[0] + dimensions[0] - boundingSizeX;
				lowerBound[1] = roiCorner[1];
				upperBound[1] = -roiCorner[1] + dimensions[1] - boundingSizeY;
				lowerBound[2] = roiCorner[2];
				upperBound[2] = -roiCorner[2] + dimensions[2] - boundingSizeZ;

				if (roiCorner[0] > bounds[i][1] + SHRINK_SIZE * 3 || (roiCorner[0] + boundingSizeX) < bounds[i][2] - SHRINK_SIZE * 3 ||
					roiCorner[1] > bounds[i][3] + SHRINK_SIZE * 3 || (roiCorner[1] + boundingSizeY) < bounds[i][4] - SHRINK_SIZE * 3)
				{
					continue;
				}

				std::cout << "Before1 : " << count << std::endl;
				generateAndSave<CTImageType>(resampledImage, resampledLabelMap, i, lowerBound, count, false);
				/*
				std::cout << "Before2 : " << count << std::endl;
				ImagePointer croppedImage = BaseAlgorithm::CropImage<CTImageType, ImageType>(resampledImage, lowerBound, upperBound);
				LabelMapPointer croppedLabelMap = BaseAlgorithm::CropImage<LabelMapType, LabelMapType>(resampledLabelMap, lowerBound, upperBound);
				for (int add = 0; add <= 0; add += 100)
				{

					// generate memory
					LabelMapPointer finalLabelMap = BaseAlgorithm::ThresholdImage<LabelMapType>(croppedLabelMap, i * 100 + 100 + add, i * 100 + 100 + add, 1);
					LabelMapPointer finalMemory = BaseAlgorithm::ThresholdImage<LabelMapType>(croppedLabelMap, 1, i * 100 + 99 + add, 1, 0);
					ImagePointer finalDistance = distanceTransform(finalLabelMap);
					std::cout << "Original phase3 : " << count << " add : " << add << std::endl;
					writeResult(croppedImage, finalMemory, finalLabelMap, finalDistance, count);

					// gauss noise or gauss smoothing
					int useSmoothing = Random(2);
					ImagePointer enhanced;
					if (useSmoothing == 1)
					{
						enhanced = gaussSmoothing(croppedImage);
					}
					else
					{
						enhanced = gaussNoise(croppedImage);
					}
					std::cout << "enhanced phase3: " << count << " add : " << add << std::endl;
					writeResult(enhanced, finalMemory, finalLabelMap, finalDistance, count);

					// flip image
					int flipAxis = Random(2);
					ImagePointer flippedImage = BaseAlgorithm::FlipImage<ImageType>(croppedImage, flipAxis);
					LabelMapPointer flippedLabelMap = BaseAlgorithm::FlipImage<LabelMapType>(finalLabelMap, flipAxis);
					LabelMapPointer flippedMemory = BaseAlgorithm::FlipImage<LabelMapType>(finalMemory, flipAxis);
					finalDistance = BaseAlgorithm::FlipImage<ImageType>(finalDistance, flipAxis);
					std::cout << "flip phase3: " << count << " add : " << add << std::endl;
					writeResult(flippedImage, flippedMemory, flippedLabelMap, finalDistance, count);

					int useRotate = Random(2);
					if (useRotate == 1)
					{
						std::cout << "rotate phase3: " << count << " add : " << add << std::endl;
						double angleX = double(Random(6) - 6 / 2) / 57.3;
						double angleY = double(Random(12) - 12 / 2) / 57.3;
						double angleZ = double(Random(16) - 16 / 2) / 57.3;
						ImagePointer rotatedImage = transformImage<ImageType>(croppedImage, angleX, angleY, angleZ);
						finalLabelMap = transformImage<LabelMapType>(finalLabelMap, angleX, angleY, angleZ, 1);
						finalMemory = transformImage<LabelMapType>(finalMemory, angleX, angleY, angleZ, 1);
						finalDistance = distanceTransform(finalLabelMap);
						writeResult(rotatedImage, finalMemory, finalLabelMap, finalDistance, count);

					}
				}
				*/
				sampled++;
			}
		} // 遍历label
	} // 遍历图像

	system("pause");
};